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Abstract 


This  paper  describes  a  method  of  obtaining  results  from  the 
simulation  of  a  finite  state  positive  recurrent  aperiodic  Markov 
chain  at  a  cost  considerably  below  the  cost  required  to  achieve  the  same 
accuracy  with  pure  random  sampling.  By  reorganizing  k  independent 
epochs  or  tours  simulated  serially  into  k  replications  simulated 
in  parallel,  one  can  induce  selected  joint  distributions  across  repli¬ 
cations  that  produce  the  cost-saving  benefits.  The  joint  distributions 
follow  from  the  use  of  rotation  sampling,  a  special  case  of  the  anti¬ 
thetic  variate  method. 

For  a  finite  state  nearest  neighbor  chain  the  paper  shows  that  even 
for  independent  parallel  replications  the  cost  of  achieving  a  specified 
accuracy  with  serial  simulation  relative  to  the  cost  for  parallel  simu¬ 
lation  has  a  lower  bound  O(k^)  as  k-^<».  When  rotation  sampling  is 
used  this  bound  is  0(k^/(ln  k)*  ).  This  lower  bound  also  holds  for 
the  more  general  finite  state  chain.  A  simulation  of  the  M/M/1 
queueing  model  with  finite  capacity  n  is  used  to  illustrate  the 
effectiveness  of  the  technique  for  selected  values  of  k,n  and 
activity  level  p  . 


Introduction 


One  ostensible  attraction  that  computer  based  simulation  offers 
to  an  experimenter  is  the  ability  to  control  the  sources  of  variation 
in  a  problem  under  study.  In  their  monograph  on  Monte  Carlo  experiments, 
Hammersley  and  Handscomb(1964)  enumerate  ways  in  which  one  can  exploit 
this  control  to  achieve  a  result  with  a  specified  accuracy  at  less  cost 
than  pure  random  sampling  would  require.  However,  that  account  prin¬ 
cipally  addresses  the  evaluation  of  univariate  and  multivariate  integrals, 
a  considerably  more  well  defined  problem  than  one  encounters  in  discrete 
event  simulation.  Although  studies  using  one  or  another  of  these  variance 
reducing  techniques  are  common  in  the  literature  on  discrete  event  simu¬ 
lation,  only  recently  has  serious  attention  been  devoted  to  developing 
methods  to  exploit  this  control  of  variation  in  conjunction  with  the 
general  classes  of  structures  encountered  there.  These  recent  develop¬ 
ments  include  elimination  of  the  list  of  scheduled  events  in  a  simulation 
containing  an  imbedded  Markov  chain  in  Hordijk,  Iglehart  and  Schassberger 
(1976),  random  number  stream  manipulation  in  multi  factor  experiments  in 
Schruben  and  Margolin  (1978),  a  perspective  on  control  variates  in 
Lavenberg  and  Welch  (1980)  and  a  reexamination  of  antithetic  variates  in 
Fishman  and  Huang  (1980). 

To  continue  this  theoretical  development,  the  present  paper  shows 
how  one  can  employ  a  special  case  of  the  antithetic  variate  technique 
called  rotation  sampling  to  exploit  the  structure  of  finite  state  Markov 
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chains  to  produce  substantial  cost  savings  in  simulation  studies  designed 
to  achieve  a  specified  level  of  statistical  accuracy.  The  generality  of 
the  approach  arises  from  the  observation  that  many  discrete  event  simuU- 
tions  have  an  underlying  Markov  chain  structure  or  a  structure  close  to 
that  of  a  Markov  chain.  Regenerative  simulation,  which  includes  Markov 
chains  as  a  special  case,  offers  one  example.  See  Crane  and  Iglehart 
(1975)  and  Fishman  (1978). 

The  proposed  technique  derives  its  cost-saving  potential  from  viewing 
the  simulation  of  k  tours  in  series  of  a  finite  (n+1)  state  positive 
recurrent  aperiodic  Markov  chain  as  equivalent  to  the  simulation  of  k 
replications  of  the  Markov  chain  in  parallel.  Although  the  marginal 
distributions  that  arise  with  the  two  alternative  formulations  are 
necessarily  the  same  for  corresponding  variables,  the  parallel  formulation 
allows  one  to  induce  joint  distributions  across  replications  that  lead  to 
a  significant  cost  saving.  The  induced  joint  distributions  follow  from 
the  use  of  rotation  sampling,  as  described  in  detail  in  Fishman  and 
Huang  (1980).  The  cost  saving  arises  in  two  ways.  Firstly,  for  fixed 
n  run  time  in  the  correlated  case  is  0(ln  k)  in  contrast  to  0(k) 
for  the  serial  simulation.  Secondly,  for  fixed  n  the  variance  of  an 
estimator  has  an  upper  bound  0(  (In  k/k)‘)  for  the  correlated  case 
compared  to  0(l/k)  for  the  serial  case.  Moreover,  the  technique 
appears  to  be  well  suited  for  use  with  the  aforementioned  proposal  of 
Hordijk,  Iglehart  and  Schassberger  (1976)  for  eliminating  the  list  of 
scheduled  events  when  simulating  a  system  representable  as  a  Markov  chain. 
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Section  1  introduces  the  Markov  chain  notation  and  specializes 
it  to  the  case  of  the  nearest  neighbor  model.  It  also  formulates  the 
experiment  of  k  indenendent  tours  or  epochs  in  series  and  then  shows 
how  one  can  reformulate  the  simulation  equivalently  as  k  independent 
replications  in  parallel.  This  reformulation  alone  shows  that  one  can 
achieve  a  substantive  cost  saving  merely  by  using  an  efficient  method 
of  generating  samples  from  the  binomial  distribution. 

Section  2  shows  how  rotation  sampling  leads  to  a  considerably  great¬ 
er  cost  saving  and  Section  3  demonstrates  how  the  results  apply  in  prac¬ 
tice.  Section  4  considers  finite  state  chains  of  a  more  general  char¬ 
acter  than  the  nearest  neighbor  model  and  extends  the  results  of  Section 
2  appropriately.  Section  5  discusses  several  unresolved  issues  which 
include  nonconvergence  to  normality,  variance  etimation,  sequential 
estimation,  infinite  state  Markov  chains  and  transient  simulation. 
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1 .  The  Nearest  Neighbor  Model 


Consider  a  positive  recurrent  aperiodic  Markov  chain  with  finite 
state  space  S  =  (0,1,. ...n)  and  (n+1)  x  (n+1)  transition  matrix 
p^  =  II  p.jjll  .  Since  all  states  0  through  n  are  regenerative, 
the  sample  paths  of  states  realized  between  successive  entries  into  a 
given  state  are  independent  and  obey  the  same  probability  law.  More¬ 
over,  a  like  property  holds  for  functions  of  these  sample  paths.  For 
expositional  convenience,  we  take  state  0  as  the  point  of  regeneration 
and  say  that  each  entry  into  this  state  completes  an  epoch  or  tour.  Let 
A.j  denote  the  reward  received  when  a  jump  occurs  from  state  i  to 

state  j,  let  N..  denote  the  number  of  transitions  from  state  i 
1  jm 

to  state  j  between  the  (m-l)st  and  mth  successive  entries  into 
state  0  (epoch  or  tour  m)  for  m  2  1  and  assume  that  observation 
begins  with  an  entry  to  state  0.  Let  u  denote  the  mean  reward  per 
epoch.  Then 


i.j=0 


*ij  ^ijm 


(1) 


is  an  unbiased  estimator  of  u  and  var  «  1/k.  Moreover,  the  number  of 
transitions  per  epoch  T  has  mean 


n 

W  =  ‘^(Pn)  =  I  E 

i,j=0 


(2) 


For  convenience  of  exposition  we  consider  the  special  but  important 
case  of  the  nearest  neighbor  Markov  chain  where 
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with  Pqq  =  0  =  1  -  Pq^ .  Section  4  considers  the  more  general  finite 
chain. 

Note  that  in  each  epoch  independent  Bernoulli  sampling  dictates  the 
branch  taken  at  each  step  or  transition.  Let  denote  a  sequence 

of  i.i.d.  random  variables  from  U(0,1)  .  Suppose  that  the  chain  is  in 
state  j  after  transition  I  -  1  and  prior  to  transition  i  .  Then  the 
reward  at  step  I  is 


where  j'  =  j'(j)  =  max{0,j-l),  j"  =  j"(j)  =  min(n,j+l)  and 


^Ca,b)^^^  "  ^ 


a  <  X  <  b 
otherwise. 


Let  T|^  and  Sj^,  denote  the  time  of  the  kth  return  to  state  0  and 
the  number  of  states  visited  on  k  such  epochs  respectively.  Then 
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E  =  ku)  and  clearly  E  Sj^  *  E  =  0(k)  .  Moreover,  the  mean  computing 
cost  of  simulating  k  epochs  is  proportional  to  E  Sj^  as  k-^®  . 

PajuxZJLtl  SJjiuLation 


Note  that  epochs  in  the  aforementioned  model  occur  in  aevial  order. 

An  alternative,  but  entirely  equivalent,  representation  in  'parallel 
facilitates  later  development  of  our  proposed  variance  reduction  method. 

In  lieu  of  running  k  simulated  epochs  in  sequence,  suppose  one  runs 
k  >  n  independent  replications  in  parallel.  Each  replication  begins 
with  an  entry  into  state  0  and  ends  with  the  next  entry  into  state  0  . 
Figure  1  gives  an  example  for  k  =  5  and  n  =  3.  Let  denote  the  number 
of  replications  in  state  i  after  transition  t  on  all  k  replications. 
One  can  then  regard  the  k  Markov  chains  operating  in  parallel  as  a  new 


Markov  chain  with  (n+l )-dimensional 


state  vector  !(,=  . 


n 

with  y  K.j  =  k  ,  after  transition  I  and  prior  to  transition  £  +  1  . 
j=0 


For  the  new  chain  an  epoch  begins  and  ends  upon  entry  into  (k,0,...,0). 


We  define  this  epoch  to  begin  at  transition  i  ~  0.  To  establish  equiva¬ 


lence  between  the  serial  and  parallel  representations,  one  needs  to  set 


PqI  =  0  =  1  -  Pqq  in  immediately  after  £  =  1  ,  so  that  states 
l,...,n  are  trone-ient, state  0  is  absorbing  and  p^  is  an  absorbing 


chain.  Note  that  the  new  state  vector  has  n  +  1  elements  and 

new  states  exist.  See  Johnson  and  Kotz  (1977,  p.  120). 

Let  i  =  l,-..,k-,  2,=  1,2,...}  denote  a  sequence  of  i.i.d. 

random  variables  from  U{0,1)  where  determines  the  branch  taken 

by  replication  i  on  step  I  .  Let  i.  .  i. 


denote  the 
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Kjj  replications  that  are  in  state  j  after  step  £  .  Then  the  number 
of  replications  that  move  from  state  j  to  j'  on  step  £  +  1  is 


_ 

’^jj'.^+l  "  ^  ^  [0,p..')^^i.  ,£+l) 

f"-l  ^JJ  j£m’  ' 


>0 


=  0 


(4) 


One  can  now  write  the  sample  mean  reward  corresponding  to  (1)  as 


n  <» 
j=0  £=0 


+  A  - 
JJ 


K.., 

Jj' 


£+1 


[(A.., 
j=0  JJ' 


^£=o'^Jj'*  ^+1 


A . y 

JJ"  ^ 


£=0 


where 


^m=l  JJ"*^  ^£=0  ’  ^m=l  ”  ^£=0  '^JJ">  ^ 


*^001  "  ^  ’ 


Sll  "  ’ 


^11  "  *^102  ^^122  "  ^ 


'  S’j'.X.+l  £+1  ^  =  1,2,...  ;  j  =  l,...,n 


Simulation  of  k  independent  replications  of  the  original  Markov 
chain  with  transition  matrix  (3)  in  terns  of  the  new  Markov  chain  for 
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l  =  1,2,...  takes  the  form 


S,t+1  Sn  '^10,J!.+1 

"  *^21,11+1 

S-,£+l  "  '^j+l,j,<t+l  ^  '^j-l.l^  ■  '^j-l,j-2,i+l  ^  ^ . 

*^n,il+l  *^n-l,!!-  ^n-1  ,n-2  ,il+l  *^n)l  *^n,n-l,£.+l 


(5) 


with  Kq^  =  0  ,  =  k  ,  =  0  for  j  =  2,...,n  and  where 

K..,  has  the  binomial  distribution  B(K..,p../)  for  j  =  l,...,n  . 

JJ  ,li+l  JJ 

Therefore,  the  distribution  of  given  depends  on  the  n 

aforementioned  binomial  distributions.  Also,  it  is  clear  that 

var  R|^  =  var  R|^  «  1/k  . 

The  computational  complexity  of  simulating  k  parallel  replications 
depends  on  two  factors,  the  cost  of  generating  I  '»  j=l,.-.,n} 

and  on  the  absorption  time  Tj^  .  If  one  use§  Bernoulli  sampling  to  generat- 
this  vector,  then  at  each  transition  the  cost  is  0(kn).  However,  algorithms  now 
exist  for  directly  sampling  Kjj/j,  and  from  the  binomial  distribution  with 

cost  0(1)  so  that  can  be  generated  at  a  cost  0(n)  .  See  Ahrens 
and  Dieter  (1980)  and  Fishman  (1978). 

Let  S|^  denote  the  total  number  of  transient  states  l,...,n  visited 
on  all  transitions  up  to  absorption.  Then  the  mean  cost  of  simulating  k 
parallel  replications  is  proportional  to  E  if  one  adopts  efficient  binomial 

sampling.  Also,  E  Sj^  <  n  E  Tj^  .  However,  note  that,  in  practice,  the  cost  of 
visiting  a  state  is  greater  in  the  parallel  case  than  it  is  in  the  serial  case. 
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More  importantly,  we  now  demonstrate  that  merely  reorganizing  a  serial  simula¬ 
tion  into  a  parallel  one  leads  to  a  cost  saving  that  grows  with  . 

Theorem  1 .  Let  V(A,B)  denote  the  cost  of  using  an  estimator  A  to 
estimate  mean  reward  relative  to  the  cost  of  using  estimator  B,  given 
var  A  =  var  B.  Then  for  fixed  n 

(i)  ET^  s  k**  (ET^)*'  . 

(ii)  V(R^.R^)  i  O(k^)  . 

See  the  Appendix  for  proof. 

The  results  of  this  theorem  imply  that  mean  absorption  time  is  sub¬ 
stantially  shorter  for  parallel  simulation  and  that  the  relat,ve  cost  of 

k 

serial  simulation  grows  with  k^  for  fixed  n  .  The  next  section  shows 
that  yet  a  more  favorable  situation  can  be  created  for  parallel  simulation 
by  inducing  an  appropriate  joint  distribution  across  the  Bernoulli  trials 
within  each  state. 
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where  Jl.  denotes  the  number  of  transitions  on  replication  j  that 

J 

occurs  in  the  cycle  beginning  with  an  exit  from  state  1  and  ending 
upon  entry  into  state  0  .  Now  observe  that  while  elements  within  a 
row  need  to  be  independent,  elements  within  a  column  need  not  be.  All 
that  is  required  within  a  column  is  that  elements  have  the  marginal  dis¬ 
tribution  Li(0,l)  . 

Consider  the  transition  1=1.  In  the  previously  described  case, 
are  independent,  so  that  var  ~  ^ar  K^22  "  *^Pio^^"PlO^  ■ 
Now  suppose  one  imposes  the  column  restriction  for  i  =  1 . k 

^•1  ■  ^11  ^ 

'  ^11  ^  Os  <  1  -  ^ 

=  if  (7) 

so  that  U.| .J , . . .  ,U|^1  each  have  U{0,1)  but  are  dependent.  This 
sampling  plan  is  a  special  case  of  rotation  sampling  and  is  discussed 
in  detail  in  Fishman  and  Huang  (1980). 

For  the  proposed  sampling  plan,  let  replication  i  jump  to  state 
0  if  0^.^  <  p^Q  and  to  state  2  if  >  P^q  .  Let  and 

denote  numbers  of  replications  that  move  from  j  to  j'  and 

j",  respectively,  on  transition  £+1  with  this  new  sampling  plan. 

Then  for  £+1  =  2,  E  kJq2  =  E  K^q2  ^  *^122  ^  ^  *^122  ^  '^^^"^10^ 
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But 


^  ^11  ~  ”  var(K^22  I  ^ii  “  k)  =  (^^^Q)(l-kp^Q) 


<  var(K^Q2  *  *^11  ■  ■  var(K 


122  ‘  ni 


k) 


(8) 


where  denotes  the  number  of  replications  in  state  j  after  transi¬ 

tion  i  and  prior  to  transition  £+1  and  kp^Q  =  kp^Q  (mod  1)  .  See 
Theorem  3  of  Fishman  and  Huang  (1980).  Note  that  under  this  new  sampling 

t  I 

plan  K^q2  ^122  have  variance  <  1/4  ,  regardless  of  k  ,  an 

attractive  feature  of  the  induced  joint  distribution^  on  U,, . U, ,  . 

II  K I 

One  can  now  extend  this  correlated  sampling  plan  to  all  subsequent 
transitions  by  grouping  replications  by  state  at  each  transition.  Recall 
that  is  the  mth  replication  on  transition  £  that  is  in  state 

I  I 

J  tor  >  0  and  m  =  l,...,Kjj^  .  Then  replication  i^.^  jumps  to 

state  j^  if  U.  .  ^  Pii'  where  for  K..  >0 
j£m 


^i  .  £  -= 

JJifn  ’j.£,i 


0  +  (mod  1 ) 


(9) 


Then  one  has: 


Theorem  2.  Assume  that  one  uses  (9)  in  a  parallel  simulation  of  the 


Fishman  and  Huang  (1981)  show  that  (8)  is  the  minimum  achievable 
variance  for  the  sum  of  k  Bernoulli  random  variables,  each  with 
the  same  marginal  distribution. 


Markov  chain  with  p  as  in  (3).  Then 


I  •  # 


(i) 

w.p. 

’  -  "Jj' 

w.p. 

ii) 

?  (1  - 

k:,  .  1/4 

(iii)  var  =  0(1)  . 


See  the  Apprndix  for  the  proof. 

The  results  in  Theorem  2  are  of  considerable  importance.  In  parti¬ 


cular,  (i)  implies  that  at  most  only  one  Bernoulli  trial  from  Bc^iCK.  P../) 

J  A.  J  J 

is  needed  for  state  j.  Thertfore,  one  can  effect  the  proposed  parallel 
sampling  plan  with  a  maximum  of  n  Bernoulli  trials.  Also  note  that  for 
large  k  the  sample  path  at  the  beginning  of  the  simulation  is  virtually 
deterministic.  Property  (ii)  establishes  the  independence  of  the  condi¬ 
tional  variance  from  k!  in  contrast  to  the  independent  case  where  this 

J  ^ 

conditional  variance  is  proportional  to  K.  .  In  particular,  observe 

J  ^ 

that  var(Kjj,^^^^|Kjjj^)  =  0  when  =  0.  Property  (iii) 

0 

removes  the  conditionality  on  and  shows  that  for  every  i  the  orders 

0  0 

of  magnitude  of  var  K..^  ,  and  var  K..„  ,  are  unaffected  by  k. 

J  J  » I  J  J  » ^ 


Simulation  of  the  Markov  chain  for  I  =  1,2...  now  takes  the  form 
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K*  =  K*  +  K* 

0,£+l  10, £+  1 

9  9 

h,i*]  “  '^2l.e+i 

K'  s  K'  4  k' 

J  +  l.j.Ul  ^-l.£ 

III  II 

“^n.il+l  "  ■  *^[1-1.0-2.11+1  *  *^0)1  ’  *^0,0-1, Ul 

with  *‘'01  '  ^  ’  *^011  '  *^11  *  ’  *^jl  '  ®  “  2,...n  and 

'^+1  t+]  Bernoulli  distribution 

8€^(KjjPjj/)  Note  that  the  principal  distinction  between  (5)  and  (10)  is 
that  (F)  relies  on  binonial  sampling,  whereas  (10)  uses  Bernoulli  sampling. 
However,  the  sample  paths  generated  on  a  single  replication  of  the  chain 
(3)  follow  the  same  probability  laws  regardless  of  the  sampling  plan. 

The  sample  mean  reward  is  now 


'^j-l.j-2.t+l 


j  =  2 . n-1 


(10) 


I 

^  j=o  e=i 


A.  .„K'.  „ 

JJ"  JJ"  ^ 


(11) 


and  the  ostensible  objective  of  analysis  Is  to  show  that  the  cost  of 
independent  serial  simulation  relative  to  the  cost  of  correlated  parallel 
simulation  for  var  R.  =  var  R*  Increases  with  k,  (k,  <  k,).  To 

^1  I  £  I 

demonstrate  this  result  we  need  to  study  T',  the  time  at  which  all  k 
correlated  replications  are  absorbed  (absorption  time). 
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Theorem  3 


and 

Then 

(i) 

(il) 

(TTi) 

(iv) 

(v) 


Let 


T|^^  =  min  (t:  I  K.^  <  r) 


j  =  1 


jt 


r  =  LI/P12J  <  » 


^kr  "  ^k  ■  ^kr 


Pjj' 


=  j.  j+i 


=  0 


£  <  j  j  =  1  .  .n. 


■  probability  of  moving  from  j  to  0  in  £  steps. 


Absorption  can  occur  at  t  if  and  only  if 


K^_1  C{(k-m,m,0,...,0); 
t-1  t-£ 

Mit  ^  I  .  n  -  I.  .  P^o  ) 


m  =  1 , . . . ,r} 


'Ot  "It  "U 


i  =  l 


n  t-j  t-e  /,•)  (,■)  t 

>  I  I  M  I  (p  -  p  Jq)  k  I 

j  =  2  £=1  i  =  l  J  ^  1  = 


no 


There  exist  b?l  and  oC(0,l)  such  that 
-b  +  k  0(.*^)  «  k  -  5  b  +  k  0(p'^)  . 


T|^^  =  0(ln  k).  w.p,  1  . 


E  T^^  =  0(ln  k)  . 
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(vi)  »  0(1)  w.p.  1  . 
(vii)  E  =  0(1n  k)  . 


See  the  Appendix  for  proof. 

Result  (i }  shows  that  total  absorption  can  occur  if  and  only  if 
K^_1  is  in  one  of  states.  Also,  result  (ii)  reveals  that  the 

rate  of  entry  into  the  absorbing  state  has  a  deterministic  component 
which  is  linear  in  k  and  a  stochastic  component  whose  magnitude  is 
independent  of  k.  Together,  results  (i)  and  (ii)  enable  one  to  establish 
parts  (iii)  through  (vii).  In  particular,  note  that  the  ratio  of  mean  absorption 

t 

times  is  E(T^)/E(lj^)  =  0(k/ln  k)  ,  clearly  favoring  correlated  parallel 
simulation  (CPS)  for  equal  k's  .  However,  it  remains  to  show  how  var 
compares  to  var  and  how  S|^.  the  number  of  states  visited  in  CPS,  com¬ 

pares  to  S|^. 

Theorem  4.  Let  c  denote  the  mean  computing  cost  of  visiting  a  state  in  the 

simulation  using  rotation  sampling  relative  to  that  cost  in  a  serial  simulation. 
For  CPS  and  fixed  n 


(i)  var  r'  •  o((ln  k/k)^) 


(ii)  ^  *  OHn  k) 


(iii) 


E  S.  var  R. 

V(Rk  ,  R,;  )  =  -V-  - -/■’ 

1  ?  E  S.  var  R, 

k2  k^ 


I 

c 


^  0(k2  /  (In  k^)^)  and  k^  <  k^ 


See  the  Appendix  for  the  proof  of  part  (1)  .  Part  (ii)  follows  from  part  (v) 
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of  Theorem  3  and  part  (iii)  follows  from  the  observation  that 
E  S.  •  var  R.  =  0  (1)  . 

kl 

Part  (iii)  is  the  result  of  ultimate  interest  in  this  study.  It  shows 

that  a  lower  bound  0(k2  /  (In  exists  on  the  cost  of 

achieving  a  specified  accuracy  with  independent  serial  tours  relative 

to  the  cost  of  achieving  the  same  accuracy  (var  R.  =  var  R'  )  with 

correlated  parallel  replications.  Moreover,  one  can  easily  show  that 

2 

this  bound  looks  like  0(k2)  as  the  accuracy  requirement  increases 
(smaller  variance).  The  next  section  demonstrates  how  these 

bounds  fare  in  practice. 
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3.  An  Illustration 

This  section  describes  a  simulation  designed  to  show  how  the 
theoretical  results  of  Section  2  fare  in  practice.  Consider  a  single 
server  queueing  system  with  independent  and  identically  distributed 
exponential  interarrival  times  with  rate  A  ,  independent  and  identi¬ 
cally  distributed  exponential  service  times  with  rate  w  >A  and 
finite  capacity  n  .  Here  jobs  that  arrive  when  n  jobs  are  already 
in  the  system  merely  go  away.  Corresponding  to  this  continuous  time 
representation  one  can  view  this  sytem  as  a  nearest  neighbor  Markov 
chain  with  p..,  =  aj/(X+w)  j=l,...,n  . 

J  «J 

As  a  figure  of  merit  we  take  mean  number  of  customers  in  system 
(u)  the  analytical  expression  for  which  appears  in  Gross  and  Harris 
(1974,  p.67).  In  the  serial  model  we  estimate  this  quantity  by 


In  the  parallel  model  using  rotation  sampling  we  estimate  u 
by 
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1 

V 

1 

X+(i) 

i 

Ti^  n  1  1 

2=1  j=l  ^ 

.  (13) 

1 

k 

1  ^ 

1 

! 

- 

Note  that  (12)  and  (13)  are  ratio  estimators  in  contrast  to  (1)  and 
(11)  which  are  Ztneor estimators.  As  a  result,  (12)  and  (13)  are  biased 
estimators  of  u  and  one  additional  inquiry  that  we  make  here  concerns 
which  estimator  has  the  smaller  bias.  Also,  it  is  of  interest  to  study 
the  extent  to  which  the  bounding  results  in  Theorem  4  apply  to  ratio 
estimation. 

t  ^ 

We  are  interested  in  observing  how  ESj^/ESj^  ,  var  Uj^/var  end 
V(C|^,C|^)  vary  with  k,n  and  p  =  X/w  .  For  convenience  and  without 
loss  of  generality  we  set  o.'  =  1  .  For  each  value  of  k  given  in  Table 
1,  1000  independent  replications  were  performed  for  each  experimental 
layout. 

Table  1 

Experimental  Layout 
k  «  2'" 


n 


'S.  II 

p\J 

1 

3 

7 

0.5 

00 

• 

II 

B 

m*2,. . . 

.9 

m=3, . . .  ,10 

0.9 

m=l ,. . .  ,8 

m*2,. . . 

.9 

m=3, . . . ,10 
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This  represents  1000  replications  on  each  of  the  48  experiments 
corresponding  to  a  unique  triplet  (k.n.p) 

Table  2  shows  selected  ratios.  Scrutinizing  these  results  one 
sees  that  E  /E  ,  var  Cj^/var  y'  and  V(C|^,y')  all  favor 
correlated  sampling.  In  particular,  for  fixed  p  and  n  ,  increas¬ 
ing  k  increases  these  ratios.  For  fixed  p  and  k  ,  increasing 
n  shows  a  degradation  in  these  ratios.  This  is  consistent  with  our 
theoretical  results.  However,  the  most  interesting  observation  occurs 
for  fixed  k  and  n  .  Here  the  ratio  var  pj^/var  p*  shows  little  effect 
due  to  changing  p  from 0.05  to  0.9.  Whether  or  not  this  observation  denotes 
an  invariance  in  the  performance  O'^  our  proposed  method  with  changes  in 
p  remains  a  topic  for  future  investigation. 

As  prevously  mentioned,  use  of  the  ratio  estimators  (12)  and 
(13)  enable  one  to  study  relative  bias.  The  ratios  (Ep|^-p)/(Ep|^-p) 
in  Table  2  show  that  correlated  sampling  always  produces  a  smaller 
absolute  bias.  However,  the  empirical  results  do  not  permit  one  to 
predict  how  the  relative  bias  changes  with  k,  n  and  p  . 

Several  other  results  not  available  in  Table  2  deserve  mention. 
Firstly,  for  the  48  triplets  considered  the  bias  in  pj^  did  not 
exceed  O.lp  and  usually  was  considerably  smaller.  Secondly,  for 
fixed  p  and  n  a  plot  of  E  Tj^ /E  Tj^  versus  k/ln  k  showed  the 

behavior  dictated  by  (iii)  of  Theorem  3.  Thirdly,  a  plot  of  var  p' 
versus  (In  k/k)^  revealed  a  linear  relationship  as  k-^  «>  confirming 
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the  applicability  of  the  lower  bound  0((ln  k/k)^)  for  this  ratio  estimator 
as  well.  Fourthly,  empirical  observation  showed  that  for  fixed  P  and 
^  ^k  ~  k -►»  .  This  too  remains  a  topic  for  future  study. 

Lastly,  we  comment  about  c,  the  mean  computing  cost  associated  with 
visiting  a  state  in  the  parallel  simulation  using  rotation  sampling 
relative  to  that  mean  cost  in  the  serial  simulation.  For  n=l  ,  c  =:  3.5 
as  k  increases;  for  n=3  and  7,0®  2.5  as  k  increases.  This  is  not 
a  serious  issue  for  three  reasons.  Firstly  increasing  k  firmly 
establishes  the  superiority  of  parallel  simulation.  Secondly,  these 
estimates  of  c  come  from  simulation  experience  in  which  many  interim 
statistics,  which  are  not  normally  collected  in  a  simulation,  were  recorded 
for  evaluation  purposes.  This  additional  collection  naturally  added  to 
computing  cost.  Thirdly,  since  no  attempt  was  made  to  optimize  the  pro¬ 
gramming  code  used,  one  expects  that  more  careful  attention  to  code  would 
reduce  c. 

4.  More  General  Markov  Chain 

We  now  consider  the  more  general  case  of  an  aperiodic  positive 
recurrent  n+1  state  Markov  chain  where  transitions  to  more  than  two 
states  are  possible,  Pqq  =  0  and  where  we  select  state  0  as  the 
regeneration  point.  Then  over  k  independent  epochs  the  sample  mean 
reward  is 

,  n  n  k 

Rl  =  i  I  I  1  a.,n.. 

^  i=0  j*0  m=l  (14) 

which  again  has  var  R|^  «  i/k  . 
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To  obtain  an  equivalent  formulation  for  k  parallel,  but  indepen¬ 
dent.  replications  we  begin  our  simulation  with  k  departures  from 


state  0  .  Let  {m.  ;  r  =  l,...s.}  denote  the  ordered  sequence 


(m.  <  m.  ,)  of  the  s  states  that  a  replication  can  enter 

J  ^  J  •  I  j  •• 


from 


J 

State  j  for  j  =  l,...,n  and  define  ^  p<  ,  Then  the 

u  ni=0 

simulation  is  for  q.  r  0: 


jm 


j.O 


1.  Move  replication  i  to  state  _  if  q^  « 

O.r  il 

for  r  »  l,...,n  and  i  ■  1 . k  »»here  ~  U(0,1)  . 

2.  Modify  p^  so  that  Pqq  ■  1  and  Pqj  ■  0  for  J  ■  1 . n  . 

3.  On  step  t  *  1  where  t  i  1  and  for  J  ■  1 , . . . ,n,  if  K-,  ^  0 

J  t  t 

move  replication  i.,  to  state  m..  if  q.  s  U. 

jtm  Jr 

f+1  ~  ^or  and 

J.r 

f  1  •  •  •  < » Sj  • 

4.  Stop  at  step  T*  =  min{£:  =  k;  £  -  1 ,2, . . , } 


The  sample  mean  reward  Is 

’  I 


T*-1 


A.. 


^  j=0  ^r»l  Vo  '^"'jr 


(15) 


If  one  imposes  the  restriction  (  9  )  ,  the  simulation  remains  as  above 
in  principle  with  replacing  K^j  in  steps  1  through  4  .  The 

sample  mean  reward  is  then 


1 

IT 


I 

=0  r- 


V'' 

^t«0  ^J"’jr  • 


Moreover,  one  has: 
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Theorem  5  . 

If  the  restriction  In  ( 

9  ) 

Is  Imposed 

,  then 

(i) 

For  r  *  1 , 

....Sj 

.mj^.£*l  “ 

-  IPJ  -  1 

w.p. 

max(P,Q)  - 

Q 

-  IQJ 

•  LPJ 

w.p. 

1  -  2  max(i 

\Q) 

-  LQJ 

-  [PJ  ♦  1 

w.p. 

max(P,Q)  - 

P 

where  P  - 

K*  ,q  .  , 

J"j.r-1 

0  ^  K'.q- 
JtVjr 

.  p 

P  (mod  1 ) 

and  Q  s  Q  (mod  1) 

(1i) 

'""‘’'ij.iti 

! 

H 

u"’ 

- 

1/4  . 

(ill) 

-  0(1). 

See  the  Appendix  for  the  proof. 

As  observed  for  the  nearest  neighbor  model,  a  simulation  using  (9) 
for  large  k  has  transitions  near  the  beginning  that  are  essentially 
deterministic.  Also,  the  random  component  of  a  transition  from  one  state 
to  another  has  the  Bernoulli  distribution,  since  examination  of  part  (1) 
of  Theorem  5  shows  that  only  two  of  the  three  outcomes  are  possible. 

$  9 

With  regard  to  the  total  absorption  time  and  var  ,  it  is  not 

difficult  to  show  that  they  behave  as  in  the  nearest  neighbor  case  as  a 
function  of  k  and  n  .  However,  5^  ,  the  number  of  states  visited, 

9  I 

calls  for  special  attention.  Whereas  one  has  t  c  ^  for  the 
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«  /II  / 

nearest  neighbor  model,  it  is  here  ^  ^  ^k  ' 

This  results  from  the  observation  that  s-  computations  are  necessary 

si 

I 

to  allocate  the  K.  replications  in  state  j  on  transition  £  . 

J 

We  summarize  these  observations  in  Theorem  6. 


Theorem  6.  For  the  more  general  finite  state  Markov  chain 


(i)  ET|^  =  0(ln  k) 

(ii)  var  r'  ^  0({ln  k/k)^) 

(iii)  ES'  •:  0(ln  k) 

(iv)  V(R|^,r')  >  0(kV  (In  k)')  . 

Proofs  follow  in  a  manner  analogous  to  those  for  earlier  theorems  and,  to 
conserve  space,  are  omitted  here. 

5.  Unresolved  Issues 

The  theoretical  results  of  Sections  1  and  2  together  with  empirical 
observations  of  Section  3  unequivocally  demonstrate  the  value  of  rotation 
sampling  for  reducing  variance  in  the  simulation  of  Markov  chains.  We 
now  address  several  ancillary  problems  reasonable  '  ^'tions  to  which  are 
necessary  before  one  can  move  totally  from  theory  to  practical  implementa¬ 
tion. 


-26- 


HoMnatiXij 

When  simulating  k  independent  epochs  in  series,  an  experimenter 
ultimately  relies  on  the  central  limit  theorem  to  derive  an  acceptable 
approximating  distribution  theory  for  confidence  intervals.  For  inde¬ 
pendent  parallel  replications,  there  is  some  plausibility  in  assuming 
normality  for  at  least  some  states  when  k  is  large.  This  follows 
from  the  observation  that  the  binomial  distribution  converges  to  the 
normal  as  k-*-oc  .  For  correlated  parallel  replications  using  (9),  the 
Bernoulli  distributions  in  each  state  for  the  sum  of  all  members  indi¬ 
cates  clearly  that  no  convergence  to  normality  exists,  regardless  of 
the  magnitude  of  k  .  In  summary,  superconvevjer.ee  of  variance  is 
achieved  at  the  loss  of  convergence  to  normality. 

t'o-ttance  E-itcmatccn 

Simulating  independent  epochs  serially  also  enables  one  to  estimate 
the  variance  needed  to  assess  accuracy.  Regrettably,  at  the  current 
writing  no  comparaole  ease  of  estimation  exists  for  k  correlated 

9 

parallel  replications.  For  as  in  (11)  a  conceivable  upper  bound 

/  / 

estimate  on  var  R  is  (nAT.  /2k)^  where  A=sup  max(  I  A  . . ,  1  , 1  A  . 

^  j  =  l,...,n 

The  value  of  this  estimate  remains  to  be  studied  in  detail.  For  the  rati 
estimator  (13),  no  comparable  proposal  is  currently  available. 
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Szquzntial  Estimation 


Although  one  should  not  minimize  the  importance  of  the  loss  of 

t 

convergence  to  normality  and  the  ability  to  estimate  var  R|^  ,  at 

least  one  issue  of  importance  in  simulation  mitigates  the  seriousness 

$ 

of  the  loss.  If  one  sets  out  to  estimate,  say,  E  to  within  an 
accuracy  of  ^6,  then  a  sequential  estimation  [irocedure  inevitably 

plays  a  role.  A  principal  result  in  this  area  is  in  Chow  and  Robbins 
(1965)  and  is  described  in  Fishman  (1978,  Ch.  2)  for  the  special  cir- 

I  /  m  \ 

cumstances  of  discrete  event  simulation.  Let  R|^'  '  denote  the  mth 
independent  rriacrofeplication  of  the  new  Markov  chain  were  each  macro¬ 
replication  consists  of  k  correlated  parallel  replications  of  the 
original  chain.  Then  one  runs  m*  macroreplications  with  outcomes 


n'(l)  r'(2)  '(m*) 


,  termination  occurring  after  the  macrorepl icatio  ■ 


(m*)  dictated  by  the  stopping  rule  associated  with  the  procedure.  By 

making  k  large  one  achieves  the  benefit  of  correlated  sampling  and 

allows  an  experimenter  to  opt  for  a  very  small  6  .  Since  as  i 0 

the  relevance  of  the  Chow-Robbins  result  for  practice  grows,  one  sees  that 

in  this  context  the  accelerated  convergence  that  correlated  replications 
/ 

induces  in  var  compensates  for  the  loss  of  normal  convergence  and  the 

I 

ability  to  estimate  var  Rj^  directly. 
finite  State  Space 

The  results  of  this  paper  clearly  show  that  the  lower  bound  on  the 
relative  benefit  of  correlated  parallel  sampling  diminishes  for  fixed 
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k  and  increasing  n  .  Although  the  seriousness  of  this  weakening 
of  the  boundary  remains  a  topic  for  further  research,  one  can  draw 
solace  from  Fable  3  which  gi'^es  the  mean  number  of  states  occupied 
per  transition  for  the  simulations  of  Section  4.  Observe  the  relatively 
slow  rate  of  increase  with  regard  to  k  . 

Table  3 

Sample  Mean  Number  of  Occupied  States  per  Transition 


XJ 

8 

16 

32 

64 

128 

256 

512 

1024 

3 

1.20 

1.40 

1.67 

1.92 

2,10 

2.20 

2.30 

7 

J  _.17  _ 

1 _ 

1.53 

1.77 

2.04 

2.38 

2.76 

3.22 

0=0.9 

3 

1.46 

1.76 

2.01 

2.18 

2.28 

2.38 

2.44 

7 

1.36 

1.68 

2.18 

2.79 

3.49 

4.05 

4.47 

4.80 

Tnan&iznt  SimutcitCon 

As  described  here,  the  simulation  we  have  in  mind  aims  at  estimating 
a  steady-state  characteristic  of  a  Markov  chain.  Alternatively  one  may 
be  interested  in  the  paths  between  states  i  and  j  and  function  of 
these  paths.  Since  the  choice  of  regenerative  simulation  via  state  0 
was  merely  a  convenience  of  exposition,  it  is  not  difficult  to  see  that 
in  the  transient  case  one  can  construct  correlated  parallel  replications 
which  have  all  the  desirable  features  given  in  Sections  3  and  5. 


1 
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APPENDIX 


Proof  of  Theorem  1 .  Let 


’  pr(movinq  from  state  j  to  state  0  in 


Since  the  parallel  replications  remain  independent,  one 


t=l  e=l 


k  E  T  . 


Then 


E  T-  -  var  T^  +  E'  T^  1  k  E  T^ 


from  which  it  follows  that 


E  Tj^  •  (E  T-)’^ 


For  part  5i)note  the  specification  that  var  R.  =  var  R, 


E  T|^  =  k  one  has  for  fixed  n 


V(Rk.R,)  = 


E  S. 


E  S. 


0(k'^ ) 


■  steps)  j  • 

has 

m=l  .2 


Since 
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Proof  of  Theorem  2 


One  can  write  K ■ ■ >  as 

J  J  .  >  I 


J.i  .■  +  ' 


K  , 

J  ■ 

'  I, 


.(u. 


D.p,,o'  i,  ~  ' 


+  I  mod  i  ' ) 


('••  i; 


jj 


j.f  .1 


r Lemi'id  5.1  nf  lishman  and  Huang  (1980),  we  hjve  equivalently 


JJ 


1=1 


'[0 


/  U*m-  1  \ 
.P'  \  y'  ' 


where  r  =  X 
1  pr(K 


JJ 


y  U  ,  (mod  1)  and  p  -  p 

'j.-.r 

1 1  .  Then 


Pq  ■-  1  -  nnn(  1  ,  K'p) 


=  ii'i n ( 1  ,f'/ p)  -  maxvO,  miO' 1  ,  K, 'p  -  1  ) 


-  nia*(0,  miml,  K*p  -  1  } 

u 


P,  --  1  -  P^  -  Pj  -  mdx(0,  min.l,  x'p  -  c  '' 


Let 


;a.j 


P  1  -  '  P  »  -  r'ax(0,  niin  1,  K'p  -  t  ) 

^  t'-o  ^ 


so  that 


pr(K 


JJ 


1) 


nia  X  ( 0 ,  m  i  n  1  1  ,  K'p  -  t  '  t  0 , 1  , .  .  .  ,  K ' 


herefore, 
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t)  =  0  0  <-  t  [K'pj  -  1 

=  IK'pJ)  =  1  -  K'p  ♦  LK'pJ  =  1  -  iPp 

pr(Kjj  .£^1  '  ♦  1)  *  ^ 

which  establishes  (i)  .  Note  that  K'..<  ,  -  [l^'pj  has  the  Bernoulli 

J  J  »  * 

distribution  Bc'iIk’p)  . 


Proof  of  (ii)  follows  immediately.  Since 
-  k'p  and  Hk'^,  v  +  '  ll^'pJ’{'  -  ^P)  ♦  (  il^'pj  ♦  1)^k'p  wc  have 

+  follows  ^ro:- 

Theorem  fishr.a''  and  liiiang  (19ftO). 


We  prove  (iii)  by  induction.  Observe  that 


...I  =  f  var(K'  ^ilK*.)  ♦  p  var  k’ 
jm..  ♦!  ,im,>*l  1  J.  '  ^jm  3 


)« 


I  II 

>"  =  j  .  J 


Since 


var  (K 


(K  p.  )(1 
J''  jm'' 


K^.P  ) 
J'^  jm' 


1/4 


(A. 3) 


f  var(K'  , |k'  ,  )  s  1/4 
jm.*>l '  Jk  ' 


(A. 4) 
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vdr{K',  ,  1  K'  )  =  0 

jm.l I  jO' 


varK^o  =0 


j  =  0,1 , . . .  ,n 


V'je  *  ^  '“''Vjr  "j-Jt* 


s  [(»jr  Kj-jj)  ’  *  (var  (A. 5) 


^0  that 


«'  ■'f'”"'' Vj.*-'*'  *  Vj. 

Where  sup  ,  1  .J  -  1  ,  .  .  ,r'  . 

beij  1  nn  1  r,,j  with  /'  =  I  ,  one  can  easily  show  that  the  magnitudes  of 


(A. 6' 


var  K  var  K  ^  ,  var  K 

jm.i  jm.?  ’  jni,3 


for  j  =  0 . r 


are  (ndependent  of  k  .  Therefore,  it  follows  hv  in-iuction  that 
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Proof  of  Theorem  3 

Let  denote  the  probability  of  total  absorption  on  transition 
t.  Then 


Qt  "  P'^^'^O.t-l  ^  St  ' 


k-1 


^m=o  ^  '  '^o.t-r  P"(S.t-l  = 


Recall  that  for  given  K-. 

J  ^ 


Sj'.t*!  "  ‘ ''je  ’’oJ'-' 


w.p.  1  -  K.^  Pjj, 


=  iKj,  Pjj'J  .  1  w.p.  K.,  Pjj,  . 


SO  that 


pr  (Kg^  =  k  I  =  m)  =  (k-m)  p^g  P'"('^l,t-1  "  '  S,t-1  " 

m  =  k  -  r,...,k  -  1 


=  0  otherwise. 

9 

Therefore,  =  0  unless  K^.,-]  C  Uk-m.m.O, ...  ,0)  ,m=l ,...  .r)  , 
which  proves  (i ) . 
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Part  (ii)  follows  from  merely  writing  down  the  transition  pattern 
into  each  state  and  regrouping  terms.  To  prove  (iii)  one  notes  that 


(a)  -1  <  .  1 

(b)  =  0  j  =  l.....n 

t  f  . .  |- 

(c)  1  -  I,  ^  Pji  -  0(i-.)  0  <  bj  1  j  =  l,...,n 

(d)  Pyl  ■  Pj,!j  i=j',  j'  +  l . 


Note  (c)  is  a  standard  result  for  transitions  from  a  transient  state  to  an 
absorbing  state  (e.g.  Cox  and  Miller,  1965).  Observe  that 


t-1 


IM 


It 


'?.  =  1 


'U, 


(1  -  I 


t-; 


'i=l 


P20  ■ 


n  t-j' 

^  I  I  M 

j=2  1=] 


y  *  ( p  ^  3  ^  -  p  ^  ^ 
^Pj'o  Py'o 


t-j 


t-i 


(i) 


r ,  n  - 1. .  pj„o) 


j=i  ii=i 


i  =  l 


where  b  -  1  and  is  finite  for  fixed  n  and  all  t  >  1  .  Therefore, 

^  ( i )  I  ^  ( i ) 

-b  +  k(l  -  [  p^p  )  .  k  -  ,  b  +  k(l  -  5;  ^  p^^  ) 

so  that  for  p  =  6^ 

-b  ♦  k  O(r')  -  k  -  .  b  *  k  0(.‘)  , 


which  proves  part  (iii) 
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Observe  that  r  <  b  +  k  0(f  and  -b  +  k  0(f  1  . 

/ 

Therefore,  0(ln  k)  i.  0  (In  k)  which  establishes  (iv)  and  (v). 

# 

Now  observe  that  for  t  s  the  number  of  replications  remaining 
in  the  transient  set  is  independent  of  k  .  Therefore,  =  0(l)w.p.  1 

t 

that  ^  =  0(ln  k),  proving  (vi)  and  (vii). 

Proof  of  Theorem  4(i) 


Let 


I  (A.  K.  +  A.  K.  .„  ) 
JJ  JJ  2  JJ  Jj"r 


Frof^  (iii)  of  Theorem  2  one  has  var  Dj  =  0(1)  .  From  (iii)  of  Tner.rer  3 

/ 

there  exists  a  constant  b^  such  that  Tj^^  <  b^  In  k.  Therefore,  one  can 
/ 

represent  as 


I 


=  ^  I  Dj  +e 

1 <ksb^ In  k 


where  t  is  a  random  variable  whose  distribution  is  independent  of  k. 
Finally 

var  r'  :  0  ((b^  In  k/k)'  )  =  0  ((In  k/k)')  , 


so 


which  establ ishes  (i ) . 
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Proof  of  Theorem  5 


Let 


where  0  q  '  p  -  1  and  L)  ~  ^^(0,1)  .  Observe  that 


K-m 

pr[L(q,p)  =  ii:]  --  I  pr[L(0,q)  =  j,  L(0,p)  =  j+m]  m  ^  0 . K  .  (A. 7) 

j  =  0 

Fror;  part  (i)  in  Theorem  2  we  have 

pr[L(0,q)  =  iKqJ  ]  =1-1^ 

pr[L(0,p)  =  IKpJ  ]  •=  1  -  Kp 

pr[L(0,q)  =  lKqJ+  1]  =  Kq 
P''[L(0,p)  =  IKpJ+  1]  =  Kp 


( 
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so  that 

pr[L(0,q)  =  L  Kq  J  ,  L(0,p)  =  L  Kp  J  ]  =  1  -  max(]^,i^) 
pr[L{0,q)  =  L  Kq  J  ,  L(0,p)  =  L Kp J  +  1]  =  max(0,  1^-1^) 
pr[L(0,q)  =  L  Kq  J  +  1.  l-(0,p)  =  L  Kp  j  +  1]  =  min(i^,i^) 

pr[L(0,q)  =  I  Kq  J  +  1 ,  L(0,p)  =  L  Kp  j  ]  =  max(0,  Kq  -  1^)  . 

Substituting  into  (A. 7)  gives 

pr[L(q,p)  =  [  Kp  J  -  I  Kq  j  -  1]  =  max(0,  1^  -  1^) 

pr[L(q,p)  =  L  Kp  J  -  L  Kq  J  ]  =  1  -  max{l^,  1^)  +  min(i^, 

pr[L(q,p)  =  I  Kpj  -  L  Kq  j  +  1]  =  max(0>  1^  -  1^)  . 

Then  partUlholas  for  q  .  q  .■'’"’j.m.  ■  =  ‘‘j.n,,  ,(*1 

u.  J,r-1  jr  jr’ 

and  K  =  K  . , 

Part  (ii)  follows  immediately  from  (A. 8)  and  part  (iii)  follows 
in  a  manner  similar  to  part  (iii)  of  Theorem  2. 
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